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ABSTRACT 


The  analysis  of  transient  signals  using  classical  techniques  is  frequently  not 
satisfactory.  The  Fourier  analysis  is  based  on  stationarity  of  the  signal,  and  transients  are 
non-stationary.  A  new  technique  for  the  analysis  of  this  type  of  signal,  called  the  Wavelet 
Transform,  was  applied  to  artificial  and  real  signals.  A  brief  theoretical  comparison 
between  the  Short  Time  Fourier  Transform  and  the  Wavelet  Transform  is  introduced.  A 
multiresolution  analysis  approach  for  implementing  the  transform  was  used.  Computer 
code  for  the  Discrete  Wavelet  Transform  was  implemented.  Different  types  of  wavelets  to 
use  as  basis  functions  were  evaluated. 
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I.  INTRODUCTION 


It  is  believed  that  different  transmitters  have  unique  signatures,  that  is,  the 
transient  or  turn  on-tum  off  response  is  unique  to  a  given  transmitter  or  at  least  to  a 
particular  make.  In  this  thesis  we  will  examine  this  issue.  The  resulting  process  algorithms 
are  also  useful  in  other  transient  signal  applications.  We  will  make  use  of  a  new  analysis 
technique  called  the  Wavelet  Transform. 

Analysis  of  stationary  signals  is  a  very  well  studied  subject.  The  researcher  can 
choose  among  a  variety  of  methods  and  apply  the  best  technique.  However,  the  study  of 
non-stationary  signals,  like  transients  and  signals  that  show  a  behavior  with  sharp 
transitions,  presents  a  challenge.  The  Fourier  Transform  method  is  not  appropriate 
because  it  uses  as  basis  function  the  complex  exponential  that  extends  over  infinite  time, 
while  transients  are  of  short  duration.  In  1946,  Gabor  [Ref.  1]  proposed  a  modification  to 
the  Fourier  Transform.  A  sliding  time-window  was  introduced  to  the  original  transform 
equation  with  the  intention  to  look  into  the  signal  by  small  portions.  Assuming  that  this 
small  portion  has  a  stationary  behavior,  the  analysis  of  the  signal  by  the  Short  Time 
Fourier  Transform  (STFT)  is  appropriate.  Other  methods  such  as  the  Instantaneous  Power 
Spectrum  and  the  Wigner-Ville  Distribution  [Ref  1 1],  have  also  been  developed  but  have 
similar  limitations  as  the  STFT. 

The  purpose  of  this  thesis  is  to  investigate  the  use  of  the  Wavelet  Transform  for 
the  estimation  /  classification  of  certain  transient  signals.  The  results  may  be  adapted  to 
problems  related  to  acoustic  and  electromagnetic  transients. 
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n.  TRANSIENT  SIGNAL  ANALYSIS 

A.  THE  SHORT  TIME  FOURIER  TRANSFORM  (  STFT  ) 

1.  The  Continuous  STFT 

The  Standard  Fourier  Transform  of  a  signal  is  defined  as 

00 

X(a)^  \  x(i)e-l°-dt .  (2.1) 

-00 

or  as 


X(a)  =  (;c(t),e^).  (2,2) 

Equations  2. 1  and  2.2  show  that  the  Fourier  coefficients  are  computed  as  inner  products 
of  the  signal  with  a  complex  exponential  basis  function  of  infinite  duration.  It  decomposes 
the  signal  into  frequency  components.  If  the  signal  components  vary  sharply  in  the  time 
domain,  after  being  transformed  these  local  irregularities  are  spread  out  over  the  entire 
frequency  domain  loosing  their  localization.  This  approach  may  not  be  suitable  if  the  signal 
of  interest  is  non-stationary.  However,  even  non-stationary  signals  can  be  assumed  to  have 
a  stationary  behavior  for  a  short  time.  Based  on  this  assumption,  Gabor  [Ref  1  ]  proposed 
the  modification  of  the  Standard  Fourier  Transform.  A  window  g(t)  of  short  duration  and 
centered  at  location  r  was  inserted  in  Equation  2. 1  as  follows 

00 

X(d)  -  J  .  ( 2.3  ) 

—00 

Equation  2.3  is  called  Short  Time  Fourier  Tran^orm  (  STFT  ).  If  the  analysis  window 
g(t}  is  Gaussian  the  STFT  is  called  the  Gabor  Transifonn.  The  analysis  window  g(t)  is 
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considered  to  be  part  of  the  specification  of  the  STFT.  The  short-time  section  ofx(t)  is  the 
product  x(t)g(t-r)  and  it  becomes  clear  that  changing  the  analysis  window  will,  in  general, 
change  all  the  short-time  sections  and  consequently  the  STFT  [Ref.  3],  A  major  drawback 
of  the  STFT  is  that  g(t),  as  a  fixed  duration  window,  is  accompanied  by  a  fixed  frequency 
resolution,  allowing  only  a  fixed  time-fi-equency  resolution  [Ref  2],  This  is  known  as  the 
Heisenberg  inequality  that  states  that  for  a  given  transform  pair 

CTr(TQ>j,  (2.4) 


where  o;  and  are  the  root  mean  square  of  the  variance  of  the  signal  over  time  and 
fi-equency  respectively.  The  variances  are  given  by 


00 


J  t^\g{t)\‘^dt 


J  lG(n)Pdin 


— OQ 


(2.5) 


(2.6) 


Hence,  either  time  or  fi-equency  resolution  can  be  enhanced  at  the  expense  of  a  poor 
resolution  in  the  other  domain  since  there  is  a  lower  bound  for  the  product  given  by 
Equation  2.4.  Note  that  the  time-frequency  resolution  is  fixed  for  the  whole 
time-frequency  plane  once  a  window  has  been  selected.  Figure  1  shows  the  tiling  of  the 

time-frequency  plane  for  the  STFT,  where  each  rectangle  represents  the  basic  resolution 
cell,  given  by  a, (To. 
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Figure  1.  Tiling  of  the  Time-frequency  Plane  for  the  STFT 


2.  The  Discrete-Time  STFT 

The  discrete-time  version  of  the  STFT  can  be  defined  as 


X{k,  w)  =  2  -  m)x(m)e  ^ 


-ffkm 


(2.7) 


The  analysis  window  g(m)  is  considered  to  be  part  of  the  specification  of  the  STFT,  and 
its  selection  dictates  the  trade-off  between  frequency  and  time  resolution.  The 
discrete-time  STFT  can  be  viewed  as  a  set  of  Discrete  Fourier  Transforms  (  DFT  )  of 
short  duration.  It  can  also  be  viewed  as  the  output  of  a  filter  bank. 

Rewriting  Equation  2.7  we  have 


X{k,  At)  =  2  [xiTn)e  ^^*”]g(A7  -  aw). 


(2.8) 


To  overcome  the  time-frequency  resolution  constraint  of  the  STFT  a  new 
technique  was  devised  to  create  a  basis  function  to  allow  a  multiresolution  analysis.  This 
new  technique  is  known  as  the  Wavelet  Transform. 

B.  THE  WAVELET  TRANSFORM 

To  obtain  a  multiresolution  analysis  we  must  allow  both  resolutions  cr,  and  to 
vary  in  the  time-frequency  plane.  When  the  analysis  is  viewed  as  a  filter  bank,  the  time 
resolution  <j,  increases  with  the  central  frequency  of  the  analysis  filters,  consequently  the 
frequency  resolution  cTq  is  inversely  proportional  to  £2,  or 

(2.10) 

where  c  is  a  constant.  The  analysis  filter  bank  is  then  composed  of  bandpass  filters  with 
proportional  bandwidth  (the  constant  Q  analysis).  Figure  3  shows  a  division  of  the 
frequency  domain  for  the  STFT  and  the  Wavelet  Transform. 


7t/  8  7t/4  tU2  n  CO 


Figure  3.  Division  of  the  Frequency-domain  for  the  STFT  and  the  Wavelet  Transform 
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The  Wavelet  transform  is  founded  on  dilations  and  translations  of  a  basis  function 
known  as  the  mother  wavelet.  The  family  derived  from  this  basis  function  is  of  the  form 

=  (2.11) 

where  a  is  the  scaling  parameter  and  b  is  the  shifting  or  translation  parameter.  If  a  is  large 
then  y/Jt)  becomes  stretched  and  corresponds  to  a  low  frequency  function  and  if  a  is 
small  y/Jt)  becomes  compressed  and  corresponds  to  a  high  frequency  function.  The 
dilations  and  translations  present  in  the  wavelet  function  enable  the  narrowing  and 
widening  of  the  time-frequency  plane.  Thus,  it  is  possible  to  adjust  to  the  characteristics  of 
the  analyzed  signal  [Ref  5].  Figure  4  shows  the  tiling  of  the  time-frequency  plane  for  the 
Wavelet  Transform.  Figure  5a  shows  a  mother  wavelet  function  while  Figures  5b,  5c  and 
5d  show  examples  of  its  dilations  and  translations. 


frequency 


Figure  4.  Tiling  of  the  Time-frequency  Plane 
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Figure  5c.  a  =  2.  b  =  1  Figure  5d.  a  =  1,  b  =  1. 

Figure  5.  Haar  Wavelet  and  Part  of  its  Family 


1.  The  Continuous  Wavelet  Transform 

The  Continuous  Wavelet  Transform  (  CWT )  can  be  defined  as 


C»Tx(a,*)  =  ^  J 


(2.12) 


where  is  given  by  Equation  2.11  and  the  term  \fja  is  the  wavelet  energy 

normalization  factor.  Equation  2.12  can  be  written  as 
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CWT:,(a,  b)  =  ■jf{\\iab{t\x{t)). 

The  inverse  of  the  Wavelet  Transform  is  given  by 


(2,13) 


00  00 

=  i  I \cWT,(a,b)^Vat^t)^.  (2.14; 

-00  0 

where  is  given  by 

~  J  fi  »  ( 2.15  ) 

0 

and  W(n)  is  the  Fourier  Transform  of  Equation  2.14  holds  only  if  the  wavelet 

function  y/Jt)  satisfies  the  admissibility  condition.  To  be  an  admissible  wavelet  function, 
the  term  <  cx3,  therefore  y/^  (t)  has  zero  mean  or  its  DC  value  is  zero  ( ^^0)  =  0  ).  The 

function  y/^  (t)  has  compact  support,  hence  it  acts  as  the  impulse  response  of  a  bandpass 
filter  that  decays  very  fast  to  zero  [Ref  2]. 

2.  The  Discrete  Wavelet  Transform 

The  parameters  a  and  b  of  the  mother  wavelet  y/^(t)  are  continuous  variables 
making  the  Wavelet  Transform  very  redundant  and  impractical.  Both  parameters  can  be 
discretized  by  sampling  to  obtain  a  set  of  discrete  wavelet  functions.  However,  if  the 
sampling  rate  is  too  high,  the  transform  is  close  to  the  continuous  transform  case  and 
hence  it  might  be  redundant.  If  it  is  too  low,  it  may  not  contain  the  necessary  amount  of 
information  for  reconstructing  the  original  signal.  The  solution  to  the  problem  may  be 
addressed  by  the  use  of  frames.  The  study  of  frames  was  introduced  by  Duffin  and 
Schaeffer  [Ref  6]  and  was  also  addressed  by  Daubechies  [Ref  7],  Their  work  provided  a 
general  theory  that  sets  the  balance  for  the  extreme  cases  of  over  and  undersampling  the 
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parameters  of  the  Short  Time  Fourier  Transform  and  the  Discrete  Wavelet  Transform  and 
proved  the  existence  of  numerically  stable  inversion  procedures  for  both  transforms. 

Let  the  sampling  lattice  be 


^  =  ^0^  (2.16) 

b  =  kboaQ,  (2.17) 

where  j ,  k  €  Z .  Equation  2. 1 1  becomes 

^ab(t)  =  ( 2.18  ) 

Substituting  2.16  and  2.17  into  2.18  we  obtain 

^jk(0  =  yifia’ot  -  kbo).  (2.19) 


The  set  of  wavelet  functions  defined  by  2.19  can  be  considered  as  a  discretized  version  of 
the  continuous  case  and  we  can  represent  any  signal  or  function  x(t)  ^L\R)  as 

•^(0  =  S  2  djJc\\fjk(t),  j,k  e  Z,  ( 2.20 ) 

J  k 

where  djji  is  defined  as 
00 

^JJc  =  J  -  kbo)cit .  (2.21) 


The  wavelet  basis  function  cannot  be  considered  a  basis  for  !*(/?);  it  is  called  z  frame.  The 
frame  does  not  satisfy  Parseval's  power  theorem  and  the  wavelet  expansion  may  lead  to 
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more  than  one  solution.  Daubechies  [Ref.  7]  showed  that,  the  recovery  oix(t)  is  possible 
only  if 


J,keZ,  (2.22) 

J  ^ 

00 

where  A  and  B  are  the  frame  bounds  for  the  frame  \ijkit)  and  \\x\\^  t  J  1^(01^^*^  The  set 

—00 

defined  by  2.19  is  considered  a  frame  only  if  V|/(?)  satisfies  the  admissibility  condition. 
Then  the  frame  bounds  are  constrained  by 

,4<7T^  I  0<A^B<o:>.  (2.23) 

bologao  J 
-00 

If  \\fjk(t)  is  inadmissible  it  leads  to  a  diverging  upper  bound  (B  =  co  ).  The  inequalities 
hold  for  any  and  b^.  Some  definitions  are  needed  before  we  proceed  : 

>  If  ^  =  1  then  the  frame  is  named  a  tight  frame; 

>•  If  the  removal  of  a  function  leads  to  an  incomplete  set,  it  is  an  exact  frame; 

>  An  exact  and  tight  frame  constitutes  an  orthonormal  basis  [Ref  7]. 

For  the  set  of  functions  defined  by  2.19  to  be  considered  a  basis  for  L^{R)  it  has 
to  satisfy  the  definitions  above  and,  therefore,  the  tightest  set  is  the  orthonormal 

set. 

Let  the  set  \^jk{f)  be  a  set  of  orthonormal  wavelets  then 


r  1.  7  =  /.  *  =  ^ 

^  0,  otherwise. 


( 2.24  ) 
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Both  indices  are  orthonomial  in  time  within  any  scale  j  and  also  across  all  other  scales 
[Ref.  2], 

For  practical  purposes  =  2  and  Zjq  =  leading  to  what  is  called  a  dyadic  or 
octave  lattice  [Ref  2,  4,  7],  Figure  6  shows  the  dyadic  or  octave  sampling  lattice. 


0.5  1  1.5  2  k 


Figure  6.  Sampling  Dyadic  Lattice  for  the  DWT 
If  we  substitute  the  values  *  2  and  =  1  in  Equation  2. 19  we  have 


^jk{t)  =  ( 2.25 

Mathematical  and  practical  perceptions  of  the  wavelet  transform  can  be  better  understood 
by  using  the  concept  of  multiresolution  to  define  the  effects  of  changing  scales  [Ref  8]. 
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3.  Multiresolution  Signal  Analysis 

Multiresolution  analysis  o^L^(R)  is  defined  as  a  sequence  of  closed  subspace  Vj  of 
^  ^  that  satisfy  the  following  properties  [Ref  9]: 

>  Containment 


>■  Scaling 

■^(0  ^  ^  0^(2f)  €  V/+1  for  any  function  x  e  L^(R); 

Shift  or  translation 

x{t)  e  Vo  <=>x(t  +  k)  e  Vo  for  any  funaion  jc  €  L^(R),  k  €  Z, 

>  A  function  (p(t)  e  Vo  exists  such  that,  for  Vm  e  Z,  the  set  defined  by 

=  2-'^(p(2-'/  —  k')  forms  an  orthonormal  basis  for  Vj.  (2.26) 

The  functions  defined  by  2.26  are  called  scaling  functions,  since  they  are  scaled  versions 
of  functions  in  L  (R)  .  From  the  above  definition  for  multiresolution  analysis,  it  can  be  said 
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that  the  function  (p(0  in  L^{R)  is  the  limit  of  the  approximations  (p/(0  €  Vj  as  j 
approaches  infinity  [Ref  5]  or 

(p(0  =  lim(p;it(0. 

y— ►<» 

The  variable  y  is  called  a  scale  factor.  For  j  >  0,  (pytCO  narrows,  and  consequently  Vj  is 
wider,  bemg  able  to  represent  finer  detail.  If  j  <0,  stretches,  Vj  is  narrower  and 

the  resolution  is  coarser.  Let  (p(0  C  Fo  be  a  scaling  function;  then  its  translations 
(p(^  -  A:)  span  Fo .  Let  (p(2?  -  A:)  be  dilations  of  (p(f  -  A:)  in  Fo .  We  know  that  Fo  c  Fj 
then  (p(2?-A)  span  V\ .  Therefore,  (p(0  can  be  presented  as  a  linear  combination  of 
translations  of  (p(2?)  as 

(p(/)  =  S  ho(k)(^{2t -k),  keZ.  ( 2.27 ) 

k 

Equation  2.27  is  known  as  dilation  or  fundamental  equation  and  ho(k)  are  the  scaling 

function  coefficients.  We  will  see  later  that  ho(k)czn  be  interpreted  as  a  FIR  filter.  The 
space  L^O^)  can  be  represented  via  a  Venn  diagram  as  shown  in  Figure  7. 


Figure  7.  Multiresolution  Representation 
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Another  way  of  graphically  representing  L^(R)  is  shown  in  Figure  8. 


Figure  8.  Multiresolution  Representation  ofL^(R) 


Note  in  Figure  8  that  the  areas  denoted  by  Wq,W\,W2,WicXc.,  are  the  difference 

spaces  between  (Fi,  Fo),(F2,  ^0,(^3,  K2),(r4,  F3),  respectively.  The  idea  can  be 
extended  to  an  arbitrary  Wj.  The  spaces  Wj  are  defined  to  be  the  orthogonal  complement 

of  the  spaces  Vj,  with  respea  to  Vj^\  [Ref  5]  that  is 

Vj^\  =  Vj  ©  Wj,  ( 2.28  ) 

(2.29) 

Where  the  symbol  "  ©  "  stands  for  direct  sum.  This  means  that  each  element  of  can 
be  written  as  the  sum  of  an  element  of  Wj  and  an  element  of  Vj  .  k  subspace  Wj 
provides  the  necessary  detail  information  to  go  from  subspace  Vj  to  Vj+\ .  The  direct  sum 

of  all  the  possibly  infinite  spaces  is  equal  to  L^(R)  [  Ref  10  ]. 


16 


Let  V|/(?),  as  defined  by  Equation  2.25,  be  an  element  of  any  given  subspace  Wj. 
The  subspace  Wj  is  itself  a  subspace  of  Vj+i .  Therefore,  V}/(?)  can  be  written  as  a  linear 
combination  of  translations  of  (p(t) 


V(0  =  S  ^  1  (k)(p(2t  -k\  keZ,  (  2.30  ) 

k 

where  hi(k)  represents  the  wavelet  coejficienis  that  can  be  used  as  a  high  pass  Finite 
Impulse  Response  (FIR)  filter.  These  coeflBcients  are  related  to  the  scaling  coefficients,  as 
they  span  two  orthogonal  spaces,  they  must  also  be  orthogonal.  Note  that  both  equations 
are  derived  fi-om  a  common  term  and  differ  only  by  the  coefficients  ho{k)  and  h\{k). 

Hence,  orthogonality  can  be  achieved  by  requiring  that 

{ho(k)M{k))  =  0.  (2.31) 

Assuming  that  the  number  of  coefficients  defining  Ao(^)  and  h  i  {k)  is  even  and 


hx{k)  =  {-\fho(N-\-k), 


( 2.32  ) 


where  N  is  the  number  of  coefficients,  then  the  inner  product  of  the  coefficients  is  equal  to 


N-\  N-\ 

Z  ho{k)hx{k)^  E  hQ{k){-\Yho{N- 1  -k). 

h=0  k=0 


(2.33) 


The  expansion  of  2.33  is  equal  to  zero.  Therefore,  the  space  L^(R)  is  spanned  by  a  linear 
combination  of  functions  in  the  subspace  Vo  and  functions  in  subspaces  IFq  to  IF  .  As  a 
consequence,  any  function  x(t)  in  L^(R)  can  be  written  as 
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(2.34) 


00  00  00 

x{t)=  2)  ak(^k{t)  +  'L  E  djk\\>jk{t). 

k=-oo  j=\  /c=-oo 

Equation  2.34  is  an  expansion  in  terms  of  scaling  and  wavelet  functions.  The  coefficients 
<3  A  and  djk  are  the  Discrete  Wavelet  Transform  coefficients.  They  completely  describe  the 

original  signal  and  are  defined  as; 


aA  =  <.x(0,(p#)>=  j  x{t)(s^k{t)dt, 

-00 


(2.35) 


00 

djk  =  {x{t\  \\ijk(f))  =  J  x{t)\^jk(t)dt . 

-00 


(2.36) 


Where  (Pa(0  and  V|/y*(f)  are  real  functions.  If  the  scaling  and  the  wavelet  functions  form 
an  orthonormal  basis,  Parseval's  power  theorem  holds.  For  the  expansion  defined  by  2.34, 
Parseval's  theorem  is 


E  +  2  \djk\'^. 

/c=-oo  Jte]  k=-<c 


(2.37) 


If  the  functions  are  orthonormal  and  have  compact  support  (  finite  duration  ),  the  desired 
time  localization  in  all  scales  is  achieved  [Ref  8,  10], 

a.  The  Scaling  Function  and  the  Scaling  Coefficients 
The  coefficients  hQ(k)  in  Equation  2.27  have  to  be  chosen  carefully  in 
order  to  generate  scaling  functions  that  have  compact  support  (finite  duration).  They  must 
satisfy  some  necessary  conditions  to  ensure  this  compact  suppon  [Ref  10.  11,  12]  The 
development  of  these  conditions  are  provided  in  Appendix  A. 
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1 .  The  sum  of  the  scaling  coefficients  must  be  equal  to  two 


A^-l 

Z 


«=0 


/zo(«)  =  2. 


(2.38) 


2.  If  a  solution  to  2.27  exists,  and  J (p{t)dt  5^  0,  and  J  dt  =  1 ,  and  an 
integer  number  of  translations  of  (p(f)  are  orthogonal  to  each  other,  as  defined  by 

((P(0(Pa(0>  =  J  q)(0q)(^  -  k)dt  =  5(A:)  = 


f  l,k  =  0 
\  0,  else 


then 


N-\ 

z 


/z=0 


hQ{n)ho{n-2k)  =  2h{k). 


(2.39) 


3.  The  sum  of  the  squares  of  the  coefficients  is  equal  to  two. 


Let  k=Q,  Equation  2.39  reduces  to 


Z  Mn)P=2. 

n=0 

4.  The  individual  sums  of  the  odd  and  even  terms 


(  2.40  ) 

of  /Jo(n)  is  equal  to  one 


iSi  +1S2  =  2, 


where 

=  Z  ko(2/7), 
and 


N-] 


82  =  X  /io(2/?  +  1). 

n=0 


(2.41) 

(2.42) 


( 2.43  ) 
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b.  Scaling  and  Wavelet  Functions  Regularity 

The  N  number  of  scaling  function  coefficients  has  to  be  even.  These 


coefficients  have  to  satisfy  the  linear  constraint  determined  by  Equation  2.38  and  the 


N 


bilinear  constraints  determined  by  Equation  2.39.  This  leaves  y-  1  degrees  of  freedom 

for  choosing  them.  That  will  guarantee  the  existence  of  the  scaling  function,  and  the  set 
derived  fi-om  it  will  be  orthonormal  [Ref  10]. 

The  degrees  of  fi'eedom  in  the  scaling  function  coefficients  determine  how 
regular  the  function  will  be.  Regularity  or  smoothness  of  a  function  is  defined  by  the 
number  of  times  that  a  function  can  be  differentiated.  In  other  words,  a  regular  scaling  or 
wavelet  function  can  be  defined  as  a  function  with  compact  support  in  time  and  reasonably 
localized  in  frequency. 

c.  The  Frequency  Domain 

The  Discrete  Fourier  transform  of  Equation  2.27  is  defined  as 


//(©)  =  X  ho(k)e 

k=-co 


(2.44) 


It  turns  out  that  a  solution  to  Equation  2.27  exists  only  if  H(0)  =  2,  or  the  DC  gain  of  the 
FIR  filter  equals  two. 

If  the  integer  translates  of  Equation  2.27  are  orthogonal,  then  the 
constraint  defined  by  Equation  2.39  is  true,  if  and  only  if 


|//(C0)|^  +  |//((o  +  n)\^  =  4. 


(245  ) 


Equation  2.45  is  equivalent  to  a  Finite  Impulse  Response  (FIR)  filter  and  is  called  a 
Quadrature  Mirror  Filter  (QMF)  [Ref  10,  18].  This  allows  the  implementation  of  the 
Discrete  Wavelet  Transform  via  a  Filter  Bank  technique. 
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4.  Filter  Banks  and  the  Discrete  Wavelet  Transform 

The  multiresolution  analysis  introduced  in  the  previous  section  can  be  implemented 
by  applying  a  technique  called  Pyramid  decomposition,  or  Mallat's  algorithm  [Ref.  13], 
The  algorithm  does  not  include  the  use  of  scaling  or  wavelet  functions;  only  the  scaling 
and  the  wavelet  coeflScients  need  to  be  considered.  We  derive  the  basic  relationship 
between  the  expansion  coefficients  at  a  lower  scale  in  terms  of  those  at  a  higher  scale. 
Reproducing  the  fundamental  equation 

(p(0  =  2  /io(/:)(p(2f  -k),  keZ.  ( 2.46 ) 

k 

Next,  we  dilate  and  translate  Equation  2.46  to  obtain 

(p(2''f  -  «)  =  2  /2o(^)cp(2(2^V  -n)-k)  = 

k 

=2  ho(k)(^(2J^^  t-2n-k).  ( 2.47 ) 

k 

Let  m  =  2n  +  or  A  =  m  -  2/7,  substituting  into  Equation  2.47 

(p(2^>  -  «)  =  2  /»o(w  -  2a2)(p(2-'"^*  t-m).  ( 2.48 ) 

m 

Let  Vj  be  the  subspace  spanned  by  (p(2^/  —  k)  then 

■^(0  ^  ^  ^(0  “  S  t~~  k^  .  (  2.49  ) 

k 

Consequently  x(t)  can  be  expressed  in  terms  of  scaling  functions  and  no  wavelets.  At  one 
immediately  lower  scale,  the  wavelets  will  provide  the  necessary  detail  not  available  in 
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that  scale  space  and  for  a  lower  resolution  the  detail  has  to  be  accounted  for.  Therefore, 
the  wavelets  have  to  be  included  in  Equation  2.49,  which  is  now  rewritten  as 


•^(0  =  S  -k)  +  'Z  djk2J'^\\/(2Jt  -  k) . 


(2.50) 


Where  the  term  2^^^  is  the  normalization  for  the  various  scales,  and  the  terms  ajk  and  djk 
are  the  discrete  wavelet  transform  coeflBcients.  Note  that  x(t)  belongs  to  I^+iand 
(^(2Jt-k)  and  \\/(2^t-k)  belong  to  .  If  we  substitute  Equation  2.26  into  Equation 


2.50  we  obtain 


x(t)  =  2  ajk%k(t)  +  S  djk\\fjk(t) ,  keZ. 

k  k 


(2.51  ) 


The  discrete  wavelet  transform  coefficients  can  be  written  as 


ajk  =  J  x{t)(pjk{t)dt . 


( 2.52  ) 


Substituting  Equations  2.26  and  2.48  into  Equation  2.52  leads  to 


Ojk  =  J  hQ{m  -  2k)2^*^^(p(2J*^t- m)dt. 


Ojk  =  2  ko(m  -  2k)  J  x(t)2^^^(p(2^*^ t -  m)dt . 


(2.53) 


Note  that  the  integrand  of  Equation  2.53  is  the  discrete  wavelet  transform  coefficient  of 
the  scaley+  1 ,  which  leads  to  the  following  relation 
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(2.54) 


^jk  =  S  ho(m  -  2k)aj+xj, . 

m 

The  analogous  relationship  for  djk  is  of  the  form 

djk  —  hxQn  —  2k')cij+ijc .  (2.55) 

m 

Therefore,  the  scaling  and  wavelet  coeflBcients  at  a  given  scale  j  are  obtained  by 
convolving  the  coeflBcients  at  that  scale  with  Ao(-/w)  and  hxi-m),  and  then  decimating  to 

produce  the  expansion  coeflBcients  at  scale  j-\.  This  is  equivalent  of  filtering  the  f**  scale 
coeflBcients  with  two  FIR  filters.  The  impulse  responses  are  hai-m)  and  hxi-m).  After 
decimation  or  down  sampling,  the  next  coarser  scaling  coeflBcients  are  obtained  [Ref  10], 

Figure  9  shows  the  implementation  of  Equations  2.52  and  2.53  for  three  levels  of 
analysis.  The  notation  LP  stands  for  the  lowpass  FIR  filter  (i.e.,  the  weights  hoi-m)  )  and 

HP  stands  for  the  highpass  FIR  filter  (i.e.,  the  weights  hxi-m)  ). 


Figure  9.  Three  levels  of  Wavelet  Transform  Analysis 
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The  implementation  of  the  algorithm  was  done  using  matrices.  An  example  of  a 
DWT  using  a  four  coefficient  wavelet  is  given  as  follows. 

Defining  a  wavelet  coefficient  matrix  A  as 


h„(0)  h„(l)  ho(2)  h„(3) 
h,(0)li,(l)  h,(2)h,(3) 

M0)h„(l)  h„(2)h<,(3) 
h,(0)h,(l)  h,(2)h,(3) 


ho(2) h„(3) 
h,(0)h,(l) 


h„(0)h„(l)  ho(2)h„(3) 
h|(0)h,(l)h|(2)h,(3) 
ho(0)h„(l) 
h,(2)h,(3) 


The  matrix  ^  is  a  sparse  matrix  composed  of  scaling  coefficients  ho(m),  and  wavelet 
coeflBcients  h  i  (m)  as  non-zero  entries.  The  blanks  indicate  zero  entries.  If  we  multiply  a 
given  signal  by  the  matrix  A,  the  first  row  generates  one  component  of  the  signal 
convolved  with  the  scaling  coefficients.  The  following  odd  numbered  rows  perform  the 
same  operation.  The  even  numbered  rows  generate  the  convolution  of  the  signal  with  the 
wavelet  coefficients.  Note  that  there  is  an  offset  by  two  every  two  rows.  This  offset  by 
two  is  the  equivalent  of  a  decimation  operation.  The  last  two  rows  wrap  around  like 
convolutions  with  periodic  boundary  conditions  [Ref  15]. 

The  DWT  consists  of  applying  the  matrix  A  repeatedly  to  a  signal  with 
starting  length  N  until  a  stopping  criterion  is  reached.  In  this  thesis,  when  the  number  of 
scaling  coefficients  was  equal  to  one,  we  stopped. 

As  an  illustration,  let  x(0  be  an  arbitrary  signal  of  length  16.  Applying  a 
DWT  to  it,  we  obtain; 
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Figure  10.  DWT  of  an  Arbitrary  Signal  x(0 

Note  that  at  each  level  matrix^  is  applied  to  half  of  the  values  of  the  previous  iterations. 
Once  the  d  coefficients  are  generated,  they  just  propagate  through  all  subsequent  levels. 
The  DWT  can  operate  only  at  signals  with  length  equal  to  a  power  of  two  (2",/?  e  Z). 

These  matrices  operations  are  generally  carried  until  there  is  only  one  DWT  scaling 
coefficient  a, , . 

The  Matlab  code  used  for  implementing  Mallat's  algoritmh  in  this  thesis  plots  the 
scales  of  the  DWT.  with  a  number  of  points  equal  to  half  of  the  total  length  of  the  signal. 
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III.  WAVELET  SIGNAL  PROCESSING 

A.  INTRODUCTION 

The  Discrete  Wavelet  Transform  was  implemented  using  the  Mallat's  Pyramid 
algorithm  described  in  Chapter  II.4.  All  programs  for  the  DWT  were  written  in  Matlab® 
[Ref  17]  and  are  given  in  Appendix  B.  The  processing  phase  was  divided  in  two  parts.  A 
simulation  phase  where  the  DWT  was  applied  to  four  artificial  communication  signals, 
and  a  testing  phase  where  the  DWT  was  applied  to  three  types  of  real  signals.  The  real 
signals  were  composed  of  transient  signals  created  by  metallic  objects,  radar  pulses,  and 
tum-on  /  tum-oflf  transients  of  push-to-talk  radios.  Six  types  of  mother  wavelet  functions 
were  used  to  process  the  signals. 

B.  MOTHER  WAVELETS 

The  Haar  wavelet  is  the  simplest  type  of  wavelet  basis  function.  It  is  very  well 
localized  in  time,  but  the  fi-equency  localization  is  poor  due  to  the  discontinuities  in  the 
time  domain.  Figure  1 1  shows  the  Haar  wavelet. 


1  - - - - - 1 

0  - - - 

-1  -  I - — _ 

- 1 - 1 - 1 _ 

0  0.5  1  1.5  2 

Figure  1 1 .  The  Haar  Wavelet 

Figure  12  shows  the  one-sided  spectrum  of  the  Haar  wavelet  and  its  dilations. 
Note  that  the  spearal  sidelobes  produce  a  poor  frequency  localization. 
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Figure  12.  One  Sided  Spectrum  of  the  Haar  Wavelet  and  its  Dilations 


The  Daubechies'  scaling  and  wavelet  functions  are  the  most  popular  functions  for 
use  in  DWT.  Daubechies'  scaling  coefficients  satisfy  all  necessary  conditions.  The  wavelets 
are  highly  regular.  They  have  compact  support  in  time  and  their  use  in  the  DWT  leads  to  a 
reasonable  localization  in  frequency.  The  regularity  of  this  type  of  wavelets  increases 
linearly  with  the  length  of  the  FIR  filter  [Ref.  2].  Daubechies'  wavelets  are  identified  by  the 
number  of  their  coefficients,  e.g.,  a  four  coefficient  wavelet  is  called  Fourth  Order 
Doubechies  or  Daubechies  four .  In  this  thesis,  the  short  notation  Daub(n)  is  used  to 
denote  the  n*  order  Daubechies  filter.  Figures  13,  14,  and  15  show  the  fourth,  twelfth  and 
twentieth  order  Daubechies'  wavelets,  respectively. 


Figure  13.  The  Fourth  Order  Daubechies  Wavelet 
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Figure  15.  The  Twentieth  Order  Daubechies  Wavelet 


Note  that  as  the  order  of  the  wavelet  is  increased  it  becomes  smoother.  However,  if 
computation  speed  is  a  factor  then  a  lower  order  Daubechies'  wavelet  might  be  preferred. 

Figures  16,  17,  and  18  show  the  frequency  response  of  Daub4,  Daub  12  and 
Daub20  wavelet  functions,  respectively.  Note  in  Figure  16  the  sidelobes  around  2.5 
radians.  Observe  that  in  the  subsequent  figures  these  sidelobes  are  gradually  becoming 
smaller,  as  the  order  of  the  wavelet  increases. 


Figure  16.  One  Sided  Spectrum  of  the  Daub4  Wavelet  and  its  Dilations 
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Figure  17.  One  Sided  Spectrum  of  the  Daubl2  Wavelet  and  its  Dilation: 
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Figure  18.  One  Sided  Spectrum  of  the  Daub20  Wavelet  and  its  Dilations 


The  Quadratic  5-spline  wavelets  are  piecewise  quadratic  functions.  These  wavelets 
are  called  the  Battle-Lemarie  bases  [Ref.  14]  and  are  referred  to  in  this  thesis  as  spline 
functions.  The  initial  scaling  function  (p(0  is  compactly  supported.  However,  the 
functions  2^^2(p(2^/-  derived  from  it  are  not  orthogonal.  The  spline  scaling  coefficients 
do  not  satisfy  all  the  constraints  described  in  Subsection  H.B.S.a.  The  use  of  the  spline 
function  in  the  DWT  leads  to  a  redundancy  in  the  scales.  Transient  signal  analysis  may, 
sometimes,  benefit  from  this  redundancy.  These  functions  have  an  arbitrarily  high  number 
of  derivatives  [Ref  7],  which  makes  them  extremely  regular.  Figure  19  shows  the  four 
coefficients  spline  function,  and  Figure  20  shows  the  one-sided  spectrum  of  its  dilations. 


Figure  19.  The  Fourth  Order  Spline  Wavelet 


Figure  20.  One  Sided  Spectrum  of  the  Spline  Wavelet  and  its  Dilations 


The  ad  hoc  wavelet  is  a  trial-and-error  wavelet  function.  The  of  purpose  of  using  it 
was  to  find  a  matching  wavelet  function.  This  function  would  be  matched  to  the  signals 
characteristics  and,  therefore,  enhances  the  detection  of  the  signal  when  wavelet 
transformed. 
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Figure  22.  One  Sided  Spectrum  of  the  ad  hoc  Wavelet  and  its  Dilations 

C  THE  SIMULATION  PHASE 

In  this  part,  four  artificial  communication  signals  created  by  Yayci  [Ref  16]  are 
used  as  test  data.  These  signals  have  sharp  transitions  in  time  that  are  exploited  via  the 
wavelet  transform.  They  are.  Binary  Phase  Shift  Keying  (BPSK),  On-Off  Keying  (OOK), 
Frequency  Shift  Keying  (FSK)  and  Quadrature  Phase  Shift  Keying  (QPSK).  They  are 
shown  in  Figures  23  through  26. 
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Figure  26.  QPSK  Signal 

1.  Application  of  the  DWT  to  Artificial  Signals 

The  DWT  was  applied  to  each  of  the  communication  signals.  The  first  experiment 
used  noise  fi-ee  signals.  Gaussian  White  Noise  (GWN)  was  added  to  the  signals  with  six 
levels  of  Signal-to-Noise-Ratio  (SNR).  The  SNR  levels  were  20  dB,  15  dB,  9  dB,  6  dB,  3 
dB,  and  0  dB.  Figures  27  to  33  present  the  decomposition  of  the  noise  free  signals.  The 
order  of  the  scales  in  these  plots  are  reversed,  i.e.,  the  first  DWT  level  appears  at  the  end. 
Only  the  last  four  scales  are  shown  because  the  other  ones  have  few  coefficients  and  little 
energy.  The  decompositions  for  the  other  levels  of  SNR  are  shown  in  Appendix  B. 

The  transition  points  of  the  signals  are  easily  seen  in  the  noise-free  case  but 
become  harder  to  identify  as  the  SNR  decreases.  Preprocessing  the  noisy  signals  using 
diffisrentiation,  Hilbert  transforms  and  a  combination  of  these  two  technicjues  did  not  seem 
to  improve  the  identification  of  the  transients.  Differentiating  a  signal  in  additive  GWN 
causes  the  scales  to  display  false  peaks.  The  increase  in  the  order  of  Daubechies'  wavelets 
improves  the  probability  of  identification  of  the  transitions.  The  ad  hoc  wavelet  perfonns 
equally  well  for  all  cases  of  GWN.  The  transitions  on  the  BPSK  and  the  OOK  signals  can 
be  identified  at  lower  SNR  relative  to  the  FSK  and  QPSK  signals.  This  is.  in  part,  due  to 
the  presence  of  sharper  transitions  in  time. 
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2.  Evaluation  of  the  DWT  for  the  Artificiai  Signals 

The  ability  of  identifying  the  transients  of  the  artificial  signals  embedded  in  noise 
was  the  performance  measure  for  the  wavelet  functions.  The  true  transient  locations  were 
known  a  priori.  Therefore,  a  simple  comparison  between  the  peaks  shown  in  the 
transform  domain  and  the  known  transition  instances  was  suflBcient  to  determine  how  well 
the  wavelet  function  performed.  The  figures  used  for  comparison  are  shown  in  Appendix 
B.  The  ad  hoc  wavelet  function  clearly  outperformed  the  other  ones.  This  wavelet  was 
chosen  by  trial-and-error,  and  it  produced  the  strongest  peaks  at  the  transition  instances.  It 
is  possible  that  its  eflBciency  will  be  degraded  when  working  with  other  types  of  signals. 
The  spline  wavelet  achieved  a  reasonable  result  down  to  6  dB  of  SNR.  The  twelfth  order 
Daubechies'  wavelet  presented  acceptable  results  for  all  signals  of  interest.  Hence,  it  was 
chosen  to  be  the  basis  function  for  the  analysis  of  the  real  signals  during  the  test  phase. 

D.  THE  TEST  PHASE 

This  section  presents  a  description  of  the  signals  used  and  their  respective 
transforms.  Only  the  last  four  scales  of  the  DWT  are  shown.  The  signals  were  divided  into 
three  groups  of  transients:  sounds  produced  by  metallic  objects,  radar  pulses,  and  tum-on  / 
turn-off  transients  produced  by  push-to-talk  radios.  These  transient  signals  are  embedded 
in  noise  with  unknown  power  levels. 

1.  Metallic  Sound  Transients 

The  first  group  of  transient  signals  was  obtained  from  the  sound  files  directory  of 
the  Digital  Signal  Processing  group  of  the  Electrical  and  Computer  Engineering 
Department  Computer  Center.  There  is  no  documentation  on  how  they  were  generated  or 
recorded.  The  length  of  the  transients  was  cut  to  the  immediate  lower  power  of  two. 
Figures  34.a  to  34.c  show  the  time  representation  of  sounds  produced  by  a  dropped 
metallic  object  on  a  metallic  surface,  a  gong,  and  a  latch  being  opened,  respectively. 
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Figure  35  shows  the  DWT  of  the  last  four  scales  of  the  metallic  sound  transients  using  a 
Daub  12  wavelet. 

Each  signal  has  a  particular  sound.  The  time  representation  of  the  gong  and  the 
dropped  metallic  object  show  that  both  signals  have  an  exponential  decaying.  The  sound 

produced  by  the  latch  is  clearly  different  from  the  other  ones,  with  no  particular  time 
behavior. 

The  DWT  of  these  signals  show  that  the  gong  and  the  dropped  metallic  object 
sounds  have  similar  characteristics  in  the  transform  domain.  This  is  due  to  their 
exponential  type  behavior  in  the  time  domain.  However,  the  latch  sound  has  a  very 
particular  transformation  and  the  transition  instances  appear  clearly  throughout  the  scales. 
A  preliminary  analysis  of  the  plots  shows  that  the  sounds  produced  by  the  dropped 
metallic  object  and  the  gong  can  be  classified  in  the  same  group  of  transients,  while  the 
sound  produced  by  the  latch  belongs  to  another  group.  It  is  possible  to  observe  these 
features  directly  from  the  time  domain  plots,  maybe  due  to  a  high  SNR.  Hence,  a  wavelet 
transform  is  not  necessary.  We  have  seen  during  the  simulation  phase  that  if  the  SNR  is 
low,  the  transition  instances  become  harder  to  identify  from  the  time  domain  plots.  Using  a 
DWT  is  a  good  technique  to  determine  these  transitions. 

Obtaining  more  metallic  sound  transients  and  transforming  them,  using  the  DWT 
for  the  purpose  of  their  classification  into  groups,  is  recommended  as  a  future  study. 

2.  Radar  Pulses 

Four  types  of  radar  pulses  were  obtained  from  two  types  of  radar  at  the  Naval 
Postgraduate  School.  The  first  radar  recorded  was  an  AN/APS  SNA  operating  at  a  center 
frequency  of  9.379  GHz,  having  a  nominal  pulse  duration  of  0.5  microseconds.  The  other 
three  types  of  radar  pulses  were  generated  by  an  AN/SPS  64  (V)  radar  operating  at  a 
center  frequency  of  9.4  GHz  and  nominal  pulse  duration  of  1.0,  0.5  and  0.15 
microseconds.  All  signals  were  heterodyned  to  500  MHz  using  a  local  oscillator  and  a 
mixer,  then  digitized  at  2  GHz  with  8  bits  per  sample.  There  are  64  pulses  per  record,  8 
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records  for  the  first  type  of  radar  pulse,  15  records  for  the  second,  8  records  for  the  third 
and  8  records  for  the  fourth.  Figures  36. a  through  36.d  show  the  first  pulse  of  the  first 
record  of  each  one  of  the  radars.  Figures  37  through  44  show  the  last  four  scales  of  the 
DWT  of  six  pulses  of  the  first  record  of  each  radar  using  the  Daub  12  wavelet. 

The  DWT  was  applied  to  the  first  six  pulses  of  each  record.  Although  the  pulses 
inside  one  record  were  generated  sequentially  by  the  same  radar,  their  DWT  scales  do  not 
present  a  particular  repeating  feature.  Therefore,  at  this  point,  we  could  not  classify  the 
transformed  radar  pulses.  Statistically,  the  number  of  radar  pulses  that  were  transformed 
was  not  enough  to  make  any  assessments.  Hence,  this  is  an  area  where  additional  research 
is  recommended. 

3.  Turn-On  /  Turn-Off  of  Push-To-Talk  Radios 

These  transient  signals  were  generated  by  the  tum-on  and  turn-off  of  six  different 
models  of  transmitters  produced  by  the  same  maker.  They  were  collected  and  recorded  by 
the  Naval  Security  Group  Activity,  Charleston,  SC.  Ten  samples  of  each  of  six 
push-to-talk  Motorola  radios  were  recorded.  The  carrier  frequency  of  all  of  them  was  set 
to  138.525  MHz.  The  signals  were  passed  through  a  1  MHz  bandwidth  filter  and  then 
digitized  with  a  sampling  frequency  of  5  MHz  and  a  center  frequency  of  1 .075  MHz. 
Figures  45.a  through  45.f  show  one  sample  of  the  tum-on  of  each  transmitter  and  Figures 
46. a  through  46.f  show  the  corresponding  transmitter  turn-off  .The  last  four  scales  of  the 
DWT  of  the  tum-on  transients,  using  a  Daub  12  wavelet,  are  shown  in  Fgures  47  and  48 
and  the  turn-off  transients  in  Figures  49  and  50. 

The  localization  of  the  exact  moment  of  the  tum-on  or  turn-off  of  the  transmitter 
can  be  extract  fi’om  the  DWT  of  these  transients.  Although  these  signals  are  sinusoidal, 
and  as  such  the  DWT  should  not  perform  well  on  them,  their  behavior  in  the  transform 
domain  is  different  before  and  after  the  switching  instant.  This  might  be  due  to  the  noise 
present  in  the  records.  Denoising  the  signals  before  applying  the  DWT  is  a  possibility  and 
deserves  further  studies.  The  transition  instances  can  also  be  observed  from  the  time 
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domain  plots  for  most  of  the  signals.  If  we  observe  Figure  45.d,  the  switching  instant  is 
difficult  to  identify,  although  if  we  look  at  Figure  48  the  transition  instance  is  observable  at 
levels  ten,  11  and  12.  A  spectrogram  of  the  record  in  question  was  performed.  Figures 
51.  a  and  51.b  show  the  result  of  the  spectrogram  using  a  hamming  window  and  without  a 
window,  respectively.  The  instant  of  the  switching  time  can  also  be  observed  in  both 
spectrograms.  The  results  of  wavelet  transforming  these  signals  are  inconclusive,  since  the 
spectrogram  produces  a  good  localization  of  the  transition  instant.  A  feature  extraction 
algorithm  based  on  the  energy  of  each  scale  for  the  purpose  of  classification  other  than 
detection  was  implemented.  The  results  were  not  satisfactory,  since  the  levels  of  energy 
for  a  given  scale  seemed  to  be  similar  for  all  records  and  all  transmitters. 
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level  13  level  12  level  11  level  10 


Figure  40.  DWT  of  Pulses  4,  5  and  6  of  Radar  AN/SPS-64  (V), 
Pulse  Duration  =1.0  microsecond 
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level  10 


level  10 


Figure  48.  DWT  of  the  First  Record  of  the  Tum-On  of  Transmitters  4  5  and  6 


Transmitter!  Transmitter  2  Transmitters 


Figure  49.  DWT  of  the  First  Record  of  the  Turn-Off  of  Transmitters  1,  2  and  3 
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Figure  51.  Spectrogram  of  Transmitter  4  First  Tum-On  Record: 

(a)  No  Windowing 

(b)  Using  a  128  Points  Hamming  Window 
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IV.  CONCLUSIONS 


A.  RESULTS 

The  identification  of  transient  signals  using  the  Discrete  Wavelet  Transform 
(DWT)  provided  mixed  results.  When  the  signals  of  interest  have  sharp  variations  in  the 
time  domain  and  the  Signal-to-Noise-Ratio  is  positive,  it  is  possible  to  identify  these 
behavior  changes  in  the  scales.  The  wavelet  used  as  a  basis  function  determines  the  overall 
performance  of  the  algorithm.  The  use  of  the  ad  hoc  wavelet  as  a  basis  function  improved 
the  identification  of  the  test  signals  in  low  SNR,  possibly  because  it  matched  the  signal 
characteristics  better  and  enhanced  them  relative  to  the  noise. 

The  DWT  of  the  sounds  produced  by  metallic  objects  indicated  that  it  might  be 
possible  to  use  the  wavelet  transform  for  classification  purposes.  Also,  the  latch  sound 
analysis  showed  that  the  exact  instant  of  the  transient  occurrences  can  be  established. 

The  question  of  what  should  be  the  scale  of  interest  was  brought  up  during  the 
simulation  and  the  test  phase.  The  higher  scales  have  more  energy  and  better  time 
resolution  and  should;  therefore,  have  more  importance  than  the  lower  scales. 

A  transient  detection  algorithm  for  the  transmitters  tum-on  /  tum-ofif  based  on  the 
energy  of  the  scales  was  implemented,  but  the  results  were  inconclusive.  The  transformed 
signals  have  the  same  total  energy  at  any  given  scale  when  compared  between  signals.  This 
was  perhaps  due  to  the  transmitters  being  manufactured  by  the  same  maker. 

The  DWT  is  simple  to  implement  and  is  very  fast.  Its  speed  depends  on  the  number 
of  scaling  function  coefiicients  that  is  being  used. 

B.  RECOMMENDATIONS  FOR  FUTURE  STUDIES 

The  implementation  of  an  automatic  transient  detection  algorithm  based  on  the 
scales  of  the  signal  is  the  goal.  Further  studies,  focused  on  finding  an  optimum  wavelet 
function,  are  necessary.  The  use  of  a  larger  number  of  actual  transient  signals  with  known 
characteristics  is  recommended  to  make  the  conclusions  statistically  more  reliable  An 
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additional  untapped  source  of  information  is  the  phase  of  the  DWT,  which  should  be 
investigated  if  the  signals  of  interest  are  complex. 
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APPENDIX  A.  DEVELOPMENT  OF  NECESSARY  CONDITIONS 


Necessary  Condition  1: 

Integrating  both  sides  of  Equation  2.27 


J  (^{f)dt  =  J  S  hQ(J^(^Q.t-k)dt 


=  'Lh„(k)\i^2t-k)dt. 


(A.1) 


Let  z  =  2t  —  k,  so  dz  =  dtl2,  substituting  into  Equation  A.  1 


J  {p{t)dt  =  5  S  hQ{k)^  (p(z)dz , 


(A.2) 


(A.3) 


The  left  hand  side  of  Equation  A.3  is  equal  to  1  and  so 


^  k 


or 


SAo(*)  =  2.  (A  4) 

k 


69 


Necessary  Condition  2: 

Substituting  Equation  2.27  in  Equation  2.39  we  obtain 


((p(0(p(?  -  ^)>  =  J  S  ho{n){^Q.t  -  hQ{m)^)(2{t  -k)-  m)dt 

n  m 

=  S  S  hQ{n)hQ{m)^ (p(2/ -  n)(^(2t -2k- m)dt  =  b{k)  .  (A.5) 

ft  ffl 


The  integral  of  Equation  2.44  is  non-zero  only  if 

2t-n  =  2t-2k-m=>n  =  2k  +  m. 

That  leads  to 

J  (p(2t  -  «)(p(2r  -2k-  m)dt  =  -m-2k).  (A.6) 

Substituting  in  equation  A.5  we  obtain 

2  ^  Ao(w)/2o(aw):“5(/7  —  m  —  2k)  =  h{k) .  (A.7) 

ft  fn  ^ 

The  Kronecker  delta  on  the  left  hand  side  exists  only  if  m=n-2k  and  Equation  A.7  is 
now 


^I,ho(n)hoin-2k)  =  m. 

^  #1 


or 


Z  ho{n)ho{n  -  2k)  =  2h{k) . 


(A.8) 
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Necessary  Condition  3: 
Reproducing  Equation  A.  8 

Zho(n)ho(n-2k)  =  28(k) 


Let  k  =  0,  and  substitute  in  Equation  A.  8 
I>!o(«)/io(n)  =  25(0)  =  2. 


(A.9) 


Necessary  Condition  4: 

Substituting  Equation  2.47  and  summing  over  k  we  get 
Z  Z  /to(«)/2o(«  -  2k)  =  2. 

k  n 


(A.  10) 


Splitting  Equation  A.  10  in  even  and  odd  terms 


Z  Z/Jo(2w)/io(2/:+2/w)  +  2/Jo(2/:+2m  +  l)/?o(2/w  + 1)  =2,  (All) 
k  L  m  m  J 


rearranging  terms 


Z  ^ho(2k-¥2m)\ho(2m)  +  Yt  Z/»o(2/:  +  2ot  + 1)  /jo(2m+ 1)  =  2  (A  12) 
»>  L  k  J  mL*  J 


Substituting  Equations  2.41  and  2.42  into  Equation  A.  12 

iSiZ  ko(2m)  +  82^  ho{2ni  +  1)  =  2. 


(A.13) 


Applying  the  same  equations  again  in  Equation  A.13 
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S\+Sl  =  2. 


(A.  14) 


Solving  for  S,  and  in  Equations  A.  10  and  A.  13  we  obtain 


*^1  =2/?o(2«)=  1 . 
n 

(A.  15) 

^2  “  2  +  1)  =  1  . 
n 

(A.  16) 
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APPENDIX  B.  DWT  OF  THE  ARTIFICIAL  SIGNALS  IN  NOISE 


haar  x10~®  daub4  xIO”^  daub12 


Figure  Bl.  DWT  of  BPSK,  SNR  =  20ciB 
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Figure  B5.  DWT  of  FSK,  SNR  =  20  dB 
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Figure  B8,  DWT  of  QPSK,  SNR  =  20  dB 
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Figure  B14.  DWT  of  FSK,  SNR  =  15  dB 
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Figure  B25.  DWT  of  BPSK,  SNR  =  6 
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6.  DWT  of  OOK,  SNR  =  3  dB 
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Figure  B40.  DWT  of  QPSK,  SNR  =  3  dB 
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APPENDIX  C.  MATLAB  CODES 


mapdn.m 

Purpose: 

Computes  a  two-dimensional  array  A  which  defines  the  mean  square  map  for  a 
dimensional  function /. 

Synopsis: 

A  =  mapdn(f,N,type). 

Description: 

Uses  a  =  wavedn(f,N,type)  to  compute  the  wavelet  transform  of/  which  is  a  sequence  of 
length  2n;  then  reorders  and  arranges  the  square  of  elements  (  l:2n  )  to  form  the 
two-dimensional  array  A,  which  has  n+1  rows  and  2(n-l)  columns;  the  volume  under  this 
surface  is  proportional  to  the  mean-square  value  of  /  using  four  mother  wavelets:  Haar, 
Daubechies,  Spline  and  the  ad  hoc. 

Algorithm: 

The  elements  of  the  wavelet  transform  a  are  recorded  according  to  the  strategy  described 
in  Ref  12,  Chapter  17,  so  that  when  the  matrix  A  is  formed  each  element  (squared)  is 
located  at  the  center  of  the  wavelet  it  represents. 

Limitations: 

Restricted  to  even  values  of  N  in  the  range  2  to  20  for  Daubechies'  wavelet  functions, 
being  N  —  2  the  call  for  the  Haar  wavelet  function  and  N  =  4  to  the  ad  hoc  wavelet 
function.  No  limitation  for  the  number  of  spline  function  coefficients. 
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Modifications: 

Addition  of  two  other  types  of  wavelet  functions  to  be  used;  spline  and  ad  hoc. 
Plotting  subroutine  appended  to  the  end  of  the  program. 

function  A  =  mapdn(f,N,type) 

% 

%  Copyright  (c)  D.  E.  Newland  1992 
%  All  rights  reserved 
% 

%  Modified  by  Jorge  Pitta,  by  permission  fi-om  the  author. 

% 

M  =  length(f); 
n  =  round(log(M)/log(2)) 
a  =  wavedn(f,N,type); 
b(l)  =  a(l);b(2)  =  a(2); 
forj  =  l;n-l 
for  k=  1:2^ 
index  =  2'j+k+N/2-l; 

while  index  >  2^(j+l), index  =  index-2'^;end 
b(index)  =  a(2^+k); 
end 
end 
a  =  b; 

forj  =  l:2'Xn-l) 

A(lJ)=a(l); 

end 

forj  =  2;n+l 
fork=  l:2^-2) 
for  m=  l;2^n-j+l) 
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A0,(k-])*2>-j+l)+m)  =  a(2^(j-2)+k); 
end 
end 
end 

A  =  A.*conj(A); 

%%  Printing  Subroutine 

subplot((n+l)/2,2,l),plot(abs(f),'r'),title(’signal-DAUB-COEFFS'), 
ylabel('fl^x)')  jXlabelCx');  grid 
hold  on 

for  h  =  2:l;n+l 

subplot((n+l)/2,2,h),plot(A(h,:),'r'),grid,ylabel(['level ' ,  num2str(h-2)]) 
end 

hold  off 
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wavedn 

Purpose: 

Wavelet  Transform. 

Synopsis: 

a  =  wavedn(f,N,type). 

Description: 

Returns  an  array  whose  elements  are  the  wavelet  transform  of  the  sequence  of  elements  in 
/  The  analysing  wavelet  has  N  even  coeflBcients. 

Algorithm: 

It  is  the  Mallat’s  Pyramid  Algorithm  adapted  to  a  circular  form. 

Limitations: 

It  is  necessary  for  the  sequence  f  to  have  2n  elements,  where  n  is  an  integer  number.  Its 
transform  a  also  has  2n  elements.  The  wavelet’s  coefficients  are  imported  from  the  file 
dcoefFs.m  (  for  the  Haar  and  Daubechies  wavelets),  dcoeffsl.m  (  for  the  Spline  wavelet ) 
and  dcoefFs2.m  ( for  the  ad  hoc  wavelet ). 

If  the  sequence/ is  complex,  then  the  transform  a  is  also  complex. 

Modifications: 

Addition  of  a  third  local  variable  called  type  that  determines  the  type  of  wavelet  function. 
Addition  of  three  IF  statements  to  use  with  the  local  variable  type, 
function  a  =  wavedn(f,N,type) 

% 

%  Copyright  (c)  D.  E.  Newland  1992 
%  All  rights  reserved 
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^  Modified  by  Jorge  Pitta,  by  permission  from  the  author. 
% 

M  =  length(f); 
n  =  round(log(M)/log(2)); 
if  type  =  =  l;c  =  dcoeffs(N);end 
if  type  =  =  2;c  =  dcoeffsl(N);end 
if  type  =  =  3;c  =  dcoefFs2(N);end 
clr  =  fliplr(c); 

for  j  =  1:2:N-1 ,  clr(j)  =  -clr(j) ;  end 
a  =  f; 

for  k  =  n:-l:l 
m  =  2^(k-l); 

X  =  [0];  y  =  [0]; 
for  i  =  1  ;m 
forj  =  1;N 
k0)  =  2*i-2+j; 

while  kO)  >  2*m ,  kQ)  =  k(j)-2*m  ;end 
end 

z  =  aCk); 

[mr,nc]  =  size(z); 
if  nc  >  1  ,  z  =  z' ;  end 
x(i)  =  c*z; 
y(i)  =  clr*z; 
end 

X  =  x/2  ;  y  =  y/2  ; 
a(l:m)  =  x; 
a(m+I;2*m)  =  y; 


dcoeffs.m 


Purpose: 

Generates  the  Haar  wavelet  coefficients  for  N  =  2  and  the  Nth  order  Daubechies’  wavelet 
coefficients  up  to  N  =  20. 

Synopsis: 
c  =  dcoeffs(N). 

Description: 

It  is  called  by  a  =  waden(f,N,type)  for  providing  the  N  Daubechies  wavelet  coefficients. 

Limitations: 

Restricted  to  even  values  of  N  in  the  range  2  to  20. 

% 

%  Copyright  (c)  D.  E.  Newland  1992 
% 

%  All  rights  reserved 
% 

function  c  =  dcoeffs(N) 

% 

ifN  =  =  2 
c  =  [l  1]; 
end 

ifN  =  =  4 

c  =  [(l+sqrt(3))/4  (3+sqrt(3))/4  (3-sqrt(3))/4  (l-sqrt(3))/4]; 
end 
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if  N  =  =  6 

q  =  sqrt(10);s  =  sqrt(5+2*q); 

c  =  [(]+q+s)/16  (5+q+3*s)/16  (5-q+s)/8  (5-q-s)/8  (5+q-3*s)/16  (l+q-s)/16]; 
end 


ifN  =  =  8 

c  =  [.3258030428051  ,1.010945715092  .892200 138246,-.039575026236,  - 
.264507 1 673 69,.0436 1 6300475,.04650360 1 07 1  ,-0. 1 4986989330] ; 
end 


ifN  =  =  10 

c  -  [.2264 1 8982583, .853943542705, 1 .024326944260,.  1 95766961 347,-.3426567 1 5382, 

-.045601 131884,.109702658642,-.008826800109,-.017791870102,.004717427938]; 

end 


ifN  =  =12 

c  =  [.157742432003,  .699503814075,  1.062263759882,  .445831322930, 
-.3 1 9986598891 ,-.  1 835 1 8064060,.  1 37888092974,.038923209708, 
-.04466374833 1  ,.00078325 1 1 52,.006756062363,-.001 523533805]; 
end 


ifN  =  =  14 

c  =  [.1 10099430746,.560791283626,1.031 148491636,.66437248221 1, 

-.2035 13822463,-.3 1683501 1281, .100846465010,.!  14003445160, 
-.053782452590;..02343994 1 565,.0 1 7749792379,.0006075 1 4996, 
-.00254790471 8,.000500226853]; 
end 
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ifN  =  =  16 

c  =  [.076955622 108,.4424672471 52,955486 150427,.8278 16532422, 
-.022385735333,-.401658632782,.000668194093,.  182076356847, 
-.024563901 046,-.06235020665 1  ,.0 1 9772 1 59296,. 0 123688448 1 9,- 
.006887719256,-.000554004548,.000955229711,-.000166137261]; 

end 


ifN=18 

c  =  [.053850349589,.344834303815,.855349064359,929545714366, 

.  1 88369549506,-.41 475 1 761802,-.  1 36953549025,.2 1 0068342279, 
.04345267546 1  ,-.095647264 120,.0003 548928 1 3,.03 1624165853, 
-.006679620227,-.006054960574,.002612967280,.000325814672, 
-.000356329759,.0000556455 1 4]; 
end 


ifN  =  =  20 

c=[.037717157593,.266122182794,.745575071487,9736281 10734,39763774177;. 

-.35333620 1 794,.2771 09878720,.  1 80127448534,.  1 3 1 602987102,.  1 0096657 1 1 96, 
.04 1 659248088,.0469698 1 4097,.005 1 00436968,.0 1 5 1 79002335,.001 973325365, 

.0028 1 7686590,-.000969947840,-.000 1 64709006,.000 1 32354366,-.00001 875841 6]; 
end 
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dcoeffsl.m 


Purpose: 

Generates  spline  wavelet  coefficients  for  any  integer  N  >  1 . 


Synopsis: 
c  =  dcoeffsl(N). 

Description: 

It  is  called  by  a  =  waden(f,N,type)  for  providing  the  N  spline  wavelet  coefficients. 

Limitations: 

None. 

% 

%  Created  by  Jorge  Pitta,  1994 

function  c  =  dcoeffsl(N) 

% 

c  =  triang(N)'; 
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dcoeffs2.ni 


Purpose: 

Generates  the  ad  hoc  wavelet  coefficients  for  any  integer  N  =  4. 


Synopsis: 

c  =  dcoeffs2(N). 

Description: 

It  is  called  by  a  =  waden(f,N,type)  for  providing  the  fourth  order  ad  hoc  wavelet 
coefficients. 

Limitations: 

Limited  to  only  a  fourth  order  ad  hoc  wavelet  function. 

% 

%  Created  by  Jorge  Pitta,  1 994 

function  c  =  dcoeffs2(N) 

% 

c  =  [.5  1.5  1]; 
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flltspec.m 


Purpose: 

Generates  the  spectrum  response  plot  of  the  Haar,  n*  order  Daubechies’,  spline  and  ad 
hoc  wavelet  function  and  their  dilations. 

Synopsis: 

filtspec. 

Description: 

It  calls  subroutines  tl.m  (  for  the  Haar  and  Daubechies*  wavelets  ),  t2.m  (  for  Spline 
wavelets ),  t3.m  ( for  the  Ad  Hoc  wavelet ). 

Created  by  Jorge  Pitta,  1994. 

F  =  menu('Wavelet','Daubechies','Spline','Ad  Hoc','Own'); 
ifF  =  =  l;tl;end 
if  F  =  =  2;t2;end 
ifF  =  =  3;t3;end 


x=[l  zeros(l,1023)];n=0:511; 

N=length(G); 

forim=l:N 

H(nn)=((-1  )^)*G(N+1  -nn); 
end 

w=filter(H,  1  pc);e=filter(G,  1  ,x); 

W=abs(ffishift(m(w,1024)));E=abs(ffishift(fft(e,1024))); 

W1=W(513:1024);E1=E(513:1024); 

H2=upsam(H,2); 
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Al=fft(H2,1024); 

Bl=fit(G,1024); 

C1=A1.*B1; 

H2G=ffishift(Cl); 

H2Gl=abs(H2G(5 13:1 024)); 

H4=upsam(H,3); 

G2=upsam(G,2); 

A2=fft(H4,1024); 

B2=fift(G2,1024); 

C2=A2.*B2.*B1; 

H4G2G=fftshift(C2); 

H4G2G 1  =abs(H4G2G(5 13:1 024)); 

H8=upsam(H,5); 

G4=upsam(G,3); 

A3=ffi(H8,1024); 

B3=f!t(G4,1024); 

C3=A3.*B3.*B2.*B1; 

H8G4G2G=fitshift(C3); 

H8G4G2G 1  =abs(H8G4G2G(5 13:1 024)); 

H16=upsam(H,7); 

G8=upsani(G,5); 

A4=fit(H16,1024); 

B4=ffi(G8,1024); 

C4=A4.*B4.*B3.*B2.*B1; 

H 1 6G8G4G2G=ffishift(C4); 

H 1 6G8G4G2G 1  =abs(H  1 6G8G4G2G(5 13:1 024)); 
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H32=upsam(H,9); 

G1 6=upsam(G,7); 

A5=ffi(H32,1024); 

B5=fift(G16,1024); 

C5=A5.*B5.*B4*B3.*B2.*B1; 

H32G 1 6G8G4G2G=fftshift(C5); 
H32G16G8G4G2Gl=abs(H32G16G8G4G2G(513:1024)); 

H64=upsam(H,l  1); 

G32=upsam(G,9); 

A6=ffi(H64,1024); 

B6=ffi(G32,1024); 

C6=A6.*B6*B5.*B4.*B3.*B2.*B1; 

H64G32G1 6G8G4G2G=fftshift(C6); 

H64G32G 1 6G8G4G2G1  =abs(H64G32G  1 6G8G4G2G(5 13:1 024)); 
m=max(H64G32G16G8G4G2Gl); 

y=Iinspace(0,pi,512);subplot(21 1) 

plot(y,H64G32Gl  6G8G4G2Gl/m;k'),hold  on 

plot(y,H32G16G8G4G2Gl/m,’k'), 

plot(y,H  1 6G8G4G2G 1  /m,'k'), 

plot(y,H8G4G2Gl/in,'k% 

plot(y,H4G2Gl/m,’k’); 

plot(y,H2Gl/m,'k’); 

plot(y,Wl/m;k’) 

pIot(y,El/m,'k'); 
hold  off 
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tl 

Purpose: 

Calls  the  m.file  dcoeffs.m  for  supplying  filtspec.m  with  the  Haar  and  Daubechies  wavelet 
coefficients. 

Synopsis: 

tl. 

Description: 

Produces  z.  pop-up  menu  on  screen  with  the  order  of  the  Daubechies'  wavelet  to  be  used, 
if  the  order  equals  two,  it  calls  Haar  wavelet  coefficients. 


Created  by  Jorge  Pitta,  1994. 


O  =  menu(’Order','2','4’,'6','8’,'10','12','14','16',’l  8','20'); 
ZZ  =  2*0; 

G  =  dcoeflFs(ZZ)/sqrt(2); 
for  ZZ  =  1  ;ZZ 

H(zz)  =  ((-l)^)*G(ZZ+l-zz); 
end 

KK  =  'Daub'; 


t2 

Purpose: 

Calls  the  m.file  dcoeffshm  for  supplying  filtspec.m  with  the  Spline  wavelet  coefficients. 

Synopsis: 

t2. 


Description: 

Produces  z  pop-up  menu  on  screen  with  the  order  of  the  Spline  wavelet  to  be  used,  if  the 
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Created  by  Jorge  Pitta,  1994.. 


O  =  menu('OrderV2V4V6','8V10V12V14V16V18','20') 
ZZ  =  2*0; 

G  =  dcoefFsl(ZZ)'/sqrt(2); 
forzz=  l.ZZ 

H(2z)  =  ((-l)^)*G(ZZ+l-zz); 
end 

KK  = 'Spline'; 

t3 


Purpose: 

Calls  the  m.£ae  dcoeffs2.m  for  supplying  filtspec.m  with  ad  hoc  wavelet  coefficients. 


Synopsis: 

t3. 


Description: 

Produces  pop-up  menu  on  screen  with  the  order  of  the  ad  hoc  wavelet  to  be  used 
However,  it  is  limited  to  an  order  4. 


Created  by  Jorge  Pitta,  1 994. 
ZZ  =  4; 

G  =  dcoeffs2(ZZ)/sqrt(2); 
for2z=  1;ZZ 

H(zz)  =  ((-l)'^)*G(ZZ+l-zz); 
end 

KK  =  -Weird'; 
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